# Файл et.top
#include "/usr/share/gromacs/top/oplsaa.ff/forcefield.itp"
[ moleculetype ]
; Name nrexcl
et 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 opls_135 1 ETH C1 1 -0.189 12.01
2 opls_135 1 ETH C2 2 -0.155 12.01
3 opls_140 1 ETH H1 3 0.0059 1.008
4 opls_140 1 ETH H2 4 0.0059 1.008
5 opls_140 1 ETH H3 5 0.0059 1.008
6 opls_140 1 ETH H4 6 0.0056 1.008
7 opls_140 1 ETH H5 7 0.0056 1.008
8 opls_140 1 ETH H6 8 0.0056 1.008
[ bonds ]
; ai aj funct b0 kb
1 2 1
1 3 1
1 4 1
1 5 1
2 6 1
2 7 1
2 8 1
[ angles ]
; ai aj ak funct phi0 kphi
;around c1
3 1 4 1
4 1 5 1
3 1 5 1
2 1 3 1
2 1 4 1
2 1 5 1
;around c2
6 2 7 1
7 2 8 1
6 2 8 1
1 2 6 1
1 2 7 1
1 2 8 1
[ dihedrals ]
; ai aj ak al funct
3 1 2 6 3
3 1 2 7 3
3 1 2 8 3
4 1 2 6 3
4 1 2 7 3
4 1 2 8 3
5 1 2 6 3
5 1 2 7 3
5 1 2 8 3
[ pairs ]
; список атомов 1-4
; ai aj funct
3 6
3 7
3 8
4 6
4 7
4 8
5 6
5 7
5 8
[ System ]
; any text here
first one
[ molecules ]
;Name count
et 1
Создание файлов для молекулярной динамики
%%bash
for i in *mdp
do
s=$(echo $i | cut -d '.' -f 1)
echo $s
grompp -f $i -c et.gro -p et.top -o et_$s.tpr
done
%%bash
for i in *tpr
do
mdrun -deffnm $i -v -nt 1
done
Следующая команда требовала дополнительный ввод поэтому вводилась вручную
%%bash
for i in *tpr
do
echo trjconv -f $i.trr -s $i -o $i.pdb
done
В PyMol были введены следующие команды для каждого из 5 pdb
delete all
load /Users/sfrv/Downloads/et_vr.tpr.pdb
mset 1-251
set cache_frames=0
mclear
viewport 210,260
mpng /Users/sfrv/Desktop/tmp/et_vr/et_vr.pdb
Из получившихся фото были сделаны гифки
from IPython.display import Image
print('Метод Андерсена')
display(Image("et_an.gif", format='png'))
print('Метод Берендсена')
display(Image("et_be.gif", format='png'))
print('Метод Нуза-Хувера')
display(Image("et_nh.gif", format='png'))
print('Метод стохастической молекулярной динамики')
display(Image("et_sd.gif", format='png'))
print('Метод "Velocity rescale"')
display(Image("et_vr.gif", format='png'))
Видно что при использовании метода Андерсена молекула практически не двигается, а при использовании метода стохастической молекулярной динамики наоборот слишком активно движется. Методы Брендсена и Нуза-Хувера слишком активно вращают молекулу вокруг своей оси. На мой взгляд наиболее приближенный к реальности результат дает метод velocity rescale.
%%bash
for i in *tpr
do
echo g_energy -f $i.edr -o ${i}_en.xvg -xvg none
done
#красный - кинетическая энергия, синий - потенциальная
import numpy as np
import matplotlib.pyplot as plt
files = ['et_an.tpr_en.xvg', 'et_be.tpr_en.xvg', 'et_nh.tpr_en.xvg', 'et_sd.tpr_en.xvg', 'et_vr.tpr_en.xvg']
methods = ['Метод Андерсена', 'Метод Берендсена', 'Метод Нуза-Хувера', 'Метод стохастической молекулярной динамики', 'Метод "Velocity rescale"']
for i in range(5):
a = np.loadtxt(f'{files[i]}')
plt.plot(a[:,0],a[:,1],'b-',a[:,0],a[:,2], 'r-')
plt.title(methods[i])
plt.xlabel('Время, пс')
plt.ylabel('Энергия, кДЖ/моль')
plt.show()
Из графиков видно, что при методе Андерсена молекула действительно обладает очень низкими энергиями. Методы velocity rescale и СМД сохраняют энергии молекулы на среднем уровне(30кДЖ/моль) и изменяют ее в достаточно широком диапазоне. Метод Брендсена сохраняет энергию на таком же уровне, однако изменяет ее в гораздо меньших пределах. Метод Нуза-Хувера поначалу изменяет энергию в больших пределах чем методы СМД и VR но постепенно остужает молекулу и разброс энергии становится сопоставим с разбросом в этих двух методах.
Рассмотрим распределение длинны связи С-С за время моделирования
%%bash
for i in *tpr
do
g_bond -f $i.trr -s $i -o bond_$i.xvg -n b.ndx -xvg none
done
#распределение длин связи в зависимости от метода
files2 = ['bond_et_an.tpr.xvg', 'bond_et_be.tpr.xvg', 'bond_et_nh.tpr.xvg', 'bond_et_sd.tpr.xvg', 'bond_et_vr.tpr.xvg']
for i in range(5):
a = np.loadtxt(files2[i])
plt.xlim((0.14, 0.165))
plt.bar(a[:,0],a[:,1], width = 0.0005)
plt.title(methods[i])
plt.xlabel('Длина связи, нм')
plt.show()
Графики показывают, что у методов с меньшими колебаниями энергии меньше и разброс длин
for i in range(5):
a = np.loadtxt(f'{files[i]}')
plt.hist(a[:,1], bins=80)
plt.title(methods[i])
plt.xlabel('Потенциальная энергия, кДж/моль')
plt.show()
Результат больше всего похожий на распределение больцмана дают методы СМД и VR. Кажется, что Метод Нуза-Хувера тоже, но это только из-за растянутой до 250 оси х.
times = [3.934, 4.062, 4.096, 4.594, 4.062]
for i in range(5):
print(f'{methods[i]}: {times[i]}')
Все методы занимают около 4 секунд на выполнение, но метод СМД занимает на 15% больше этого времени, а метод Андерсена оказался незначительно быстрее всех. На мой взгляд при малых температурах системы хорошо подойдут методы Андерсена и Нуза-Хувера, а при высоких -- методы СМД и VR(мне кажется он лучше описывает поведение молекулы, распределение потенциальных энергий наиболее приближено к Больцмановскому, да и к тому же он быстрее СМД).